function [P s coeffs]=specw2(x,dt,nw,opts);

%function [P s P95 Pc Pprefit]=specw(x,dt,nw,opts);
%
%INPUT
%
%x     equally-spaced time series
%dt    time spacing
%nw    number of windows (for mtm), nw=1 gives a periodogram
%opts: options
%      1. Prewhitening:
%         opts.p, number of autoregressive coefficients (default 1)
%         opts.q, number number of moving-average coefficients (default 0)
%      2. Postwhitening:
%         opts.smoothing; choose from 'gaussian' (default),
%         'lowess', or others available through smoothdata.m
%      3. opts.plot: 1=yes, 0=no (default).
%
%OUTPUT
%
%P     Squared fourier coefficients (nw=0) or power spectral
%      density (nw>0)
%s     Frequency
%P95   95% confidence interval as a function of s
%Pc    Background continuum estimate
%R     Residual, computed as log(P)-log(Pc).

if nargin<3,
    dt=1;
    nw=4;
    %simulate a time-series
    mdlV=1;
    mdl = arima('MA',[0.6],'AR',[0.5 0.2],'Variance',mdlV,'Constant',0);
    x=simulate(mdl,10^4,'NumPaths',1);
end;
if nargin<4,
    opts.p=1;
    opts.q=0;
    opts.smoothing='gaussian';
    opts.plot=1; 
else,
    if ~isfield(opts,'p'), opts.p=1; end;
    if ~isfield(opts,'q'), opts.q=0; end;
    if ~isfield(opts,'smoothing'), opts.smoothing='gaussian'; end;
    if ~isfield(opts,'plot'), opts.plot=0; end;    
end;

%Expected relationships based on the number of degrees of freedom
k=2*nw-1;
coeffs.dof=2*k;

alpha1=0.05*2*k/length(x);
coeffs.conf=log(gaminv(1-alpha1,k,1/k)); %95 confidence level relative to the log of the background continuum

alpha2=0.01*2*k/length(x);
coeffs.confb=log(gaminv(1-alpha2,k,1/k)); %99 confidence level for multiple testing regime (Bonferroni correction)

alpha3=1-(1-0.05)^((2*k)/length(x));
coeffs.confc=log(gaminv(1-alpha3,k,1/k)); %95 confidence level for
                                          %multiple testing regime
                                          %(Sidak correction)

coeffs.Vexp=psi(1,coeffs.dof/2);  %expected variance of log(P)-log(Pc)
coeffs.bias=psi(0,k)+log(2)-log(k*2); %expected bias in the mean estimate of log(Pc)

%spectra of raw record
x=detrend(x);
[dummy s]=pmtmPH(x,dt,nw,0); 
P.initial=log(dummy)-mean(log(dummy))+coeffs.bias; 

%prewhiten using an ARIMA model
[xpre dummy]=prewhitenPH(x,opts.p,opts.q);
coeffs.p=dummy.p;
coeffs.q=dummy.q;
coeffs.v=dummy.v;
[dummy s]=pmtmPH(xpre,dt,nw,0); 
P.prewhiten=log(dummy)-(mean(log(dummy))-coeffs.bias); 
coeffs.varP=var(P.prewhiten);

%implied prewhitened spectral fit
P.pre_fit=log(arma_spec(s,dt,coeffs.p,coeffs.q,coeffs.v));

